home *** CD-ROM | disk | FTP | other *** search
/ Dr. Windows 3 / dr win3.zip / dr win3 / PROGRAMR / GSRC208A.ZIP / GAUSSW.C < prev    next >
C/C++ Source or Header  |  1993-08-26  |  19KB  |  600 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*         GWSIM - Simulation        */
  11. /*        MS-WINDOWS front end       */
  12. /*                                   */
  13. /*          Gauss reduction          */
  14. /*            dialog box             */
  15. /*                                   */
  16. /*           QuickC/WIN 1.0          */
  17. /*                                   */
  18. /*   (include here compilers that    */
  19. /*   compiled GWSIM successfully)    */
  20. /*                                   */
  21. /*************************************/
  22.  
  23. #include <windows.h>
  24. #include <math.h>
  25. #include <string.h>
  26. #include "globals.h"
  27. #include "gep2.h"
  28. #include "strtbl.h"
  29.  
  30. #define ALMOST_ZERO 1.0e-7
  31.  
  32. GLOBALHANDLE    hRSto;                    /* handle to memory block w/ rstoi        */
  33. GLOBALHANDLE    hSto;                    /* handle to memory block w/ stoi        */
  34. GLOBALHANDLE    hMl;                    /* handle to memory block w/ ml            */
  35. GLOBALHANDLE    hLm;                    /* handle to memory block w/ lm            */
  36. GLOBALHANDLE    hLD;                    /* handle to memory block w/ ld            */
  37. GLOBALHANDLE    hMetn;                    /* handle to memory block w/ metn        */
  38. GLOBALHANDLE    hStepn;                    /* handle to memory block w/ stepn        */
  39. GLOBALHANDLE    hStepst;                /* handle to memory block w/ stepst        */
  40. GLOBALHANDLE    hImet;                    /* handle to memory block w/ imet        */
  41. GLOBALHANDLE    hRs;                    /* handle to memory block w/ rs            */
  42. GLOBALHANDLE    hLoo;                    /* handle to memory block w/ loo        */
  43. float            huge *rstoi;            /* pointer to work with rstoi array        */
  44. float             huge *ml;                /* pointer to work with ml array        */
  45. float             huge *lm;                /* pointer to work with lm array        */
  46. float             huge *ld;                /* pointer to work with ld array        */
  47. int                huge *imet;                /* pointer to work with imet vector        */
  48. int                huge *loo;                /* pointer to work with loo array        */
  49. int        (huge *rs)[MAX_STEP][MAX_MOL];        /* pointer to work with rs array        */
  50. char    (huge *metn)[NAME_L];            /* pointer to work with metname array    */
  51. char    (huge *stepn)[NAME_L];            /* pointer to work with metname array    */
  52. char    (huge *stepst)[256];            /* pointer to work with stepst array    */
  53. int ur[MAX_MET];                        /* permut. on metabolites  u -> g         */
  54. int uc[MAX_MET];                        /* permutations on steps   u -> g        */
  55.  
  56. void rowsw( int r1, int r2, int c );
  57. void colsw( int c1, int c2, int r );
  58. void rowsw_m( int r1, int r2, int c );
  59. void colsw_m( int c1, int c2, int r );
  60. int emptyrow( int rw, int nsteps);
  61. void invml11( void );
  62. void calc_ld( void );
  63. void init_moiety( void );
  64. void lindep( void );
  65. int gauss( void );
  66. void more_ext( void);
  67. void int_ord( void );
  68. void ext_ord( void );
  69. void initreds( void );
  70. void virt_step( void );
  71. int reduce( int swapnames );
  72.  
  73. #pragma alloc_text( CODE14, rowsw, colsw, rowsw_m, colsw_m, emptyrow, invml11, calc_ld, init_moiety, lindep, gauss, int_ord, ext_ord, more_ext, initreds, virt_step, reduce )
  74.  
  75. void rowsw( int r1, int r2, int c )
  76. {
  77.  register int j;
  78.  float dummy;
  79.  double dumb;
  80.  int dum;
  81.  
  82.  for(j=0;j<c;j++)
  83.  {
  84.   dum         = stoi[r1*MAX_MET + j];                      /* switch rows in stoi and */
  85.   stoi[r1*MAX_MET + j] = stoi[r2*MAX_MET + j];
  86.   stoi[r2*MAX_MET + j] = dum;
  87.   dummy        = rstoi[r1*MAX_MET + j];                    /* switch rows in rstoi &  */
  88.   rstoi[r1*MAX_MET + j] = rstoi[r2*MAX_MET + j];
  89.   rstoi[r2*MAX_MET + j] = dummy;
  90.  }
  91.  j      = ur[r1];                                 /* switch rows in ur and   */
  92.  ur[r1] = ur[r2];
  93.  ur[r2] = j;
  94. }
  95.  
  96. void colsw( int c1, int c2, int r )
  97. {
  98.  register int j;
  99.  float dummy;
  100.  int dum;
  101.  
  102.  for( j=0; j<r; j++ )
  103.  {
  104.   dum         = stoi[j*MAX_MET + c1];                      /* switch cols in stoi and */
  105.   stoi[j*MAX_MET + c1] = stoi[j*MAX_MET + c2];
  106.   stoi[j*MAX_MET + c2] = dum;
  107.   dummy        = rstoi[j*MAX_MET + c1];                    /* switch cols in rstoi &  */
  108.   rstoi[j*MAX_MET + c1] = rstoi[j*MAX_MET + c2];
  109.   rstoi[j*MAX_MET + c2] = dummy;
  110.  }
  111.  j      = uc[c1];                                 /* switch cols in uc too   */
  112.  uc[c1] = uc[c2];
  113.  uc[c2] = j;
  114. }
  115.  
  116. void rowsw_m( int r1, int r2, int c )             /* switch rows in ml       */
  117. {
  118.  register int j;
  119.  float dummy;
  120.  
  121.  for(j=0;j<c;j++)
  122.  {
  123.   dummy     = ml[r1*MAX_MET + j];
  124.   ml[r1*MAX_MET + j] = ml[r2*MAX_MET + j];
  125.   ml[r2*MAX_MET + j] = dummy;
  126.  }
  127. }
  128.  
  129. void colsw_m( int c1, int c2, int r )             /* switch cols in ml       */
  130. {
  131.  register int j;
  132.  float dummy;
  133.  
  134.  for(j=0;j<r;j++)
  135.  {
  136.   dummy     = ml[j*MAX_MET + c1];
  137.   ml[j*MAX_MET + c1] = ml[j*MAX_MET + c2];
  138.   ml[j*MAX_MET + c2] = dummy;
  139.  }
  140. }
  141.  
  142. int emptyrow( int rw, int nsteps)
  143. {
  144.  register int j, ct;                              /* ct counts entries != 0  */
  145.  
  146.  for( ct=0, j=0; j<nsteps; j++)
  147.   if ( fabs( rstoi[rw*MAX_MET + j] ) > ALMOST_ZERO ) ct++;
  148.  return ( ct );
  149. }
  150.  
  151. void invml11( void )
  152. {                                                 /* as ml is lower triangul.*/
  153.  register int i,j,k;                              /* inverting is simply to  */
  154.  float acum;                                      /* forward-substitute with */
  155.                                                   /* the identity matrix     */
  156.  for(i=0;i<nmetab; i++) ml[i*MAX_MET + i] = 1 ;
  157.  for( j=0; j<indmet; j++)
  158.   for( k=0; k<indmet; k++)
  159.   {
  160.    for( acum=0, i=0; i<k; i++ )
  161.     acum += ml[k*MAX_MET + i] * lm[i*MAX_MET + j];
  162.    lm[k*MAX_MET + j] = ( (k==j) ? (float) 1.0 : - acum ) / ml[k*MAX_MET + k];
  163.   }
  164. }
  165.  
  166. void calc_ld( void )
  167. {
  168.  register int i,j,k;
  169.  
  170.  for( i=indmet; i<nmetab; i++ )                   /* first calculate L0      */
  171.   for( j=0; j<indmet; j++ )                       /* which is                */
  172.   {
  173.    ld[i*MAX_MET + j] = 0;
  174.    for( k=0; k<indmet; k++ )                      /*                 -1      */
  175.     ld[i*MAX_MET + j] += ml[i*MAX_MET + k] * lm[k*MAX_MET + j];              /* L0 = L21 * (L11)        */
  176.   }
  177.  for( i=indmet; i<nmetab; i++ )                   /* Now map the lin. dep.   */
  178.  {
  179.   for( j=0; j<indmet; j++ )                       /* which is                */
  180.    ld[i*MAX_MET + j] = -ld[i*MAX_MET + j];                          /* -L0                     */
  181.   for( j=indmet; j<nmetab; j++ )                  /* concateneted with the   */
  182.    ld[i*MAX_MET + j] = (i == j)
  183.     ? (float) 1.0 : (float) 0.0;                  /* m0 -> m part of Id.     */
  184.  }
  185. }
  186.  
  187. void init_moiety( void )
  188. {
  189.  register int i,j;
  190.  
  191.  for( i=indmet; i<nmetab; i++)
  192.  {
  193.   moiety[i] = (float) 0.0;
  194.   for( j=0; j<nmetab; j++ )
  195.    moiety[i] += ld[i*MAX_MET + j] * xu[j];
  196.  }
  197. }
  198.  
  199. void lindep( void )
  200. {
  201.  invml11();                                       /* invert ml11 into lm     */
  202.  calc_ld();                                       /* calculate L2 * -L1'     */
  203.  init_moiety();                                   /* setup the moiety conc.  */
  204. }
  205.  
  206. int gauss( void )
  207. {
  208.  register int i, j;
  209.  int k, flag;
  210.  float m;
  211.  
  212.  for( k=0; k<nsteps+1; k++)
  213.  {
  214.   for( flag=0, i=k; i<nmetab; i++ )
  215.    if( fabs(rstoi[i*MAX_MET + k]) > ALMOST_ZERO )
  216.    {
  217.     flag = 1;
  218.     break;                                        /* suitable divisor        */
  219.    }
  220.   if ( flag )
  221.   {
  222.    if ( i != k )
  223.    {
  224.     rowsw( k, i, nsteps );                        /* switch row k with i     */
  225.     rowsw_m( k, i, nmetab );
  226.    }
  227.   }
  228.   else
  229.   {
  230.    do
  231.    {
  232.     for( j=k+1; j<nsteps; j++ )
  233.      if( fabs(rstoi[k*MAX_MET + j]) > ALMOST_ZERO )
  234.      {
  235.       flag = 1;
  236.       break;                                      /* suitable value          */
  237.      }
  238.     if ( flag && ( j != k ) )
  239.     {
  240.      colsw( j, k, nmetab );                       /* switch col j with k     */
  241.      colsw_m( j, k, nmetab );
  242.     }
  243.     if( !flag )
  244.     {
  245.      for( i=0; i<nmetab; i++ )                    /* get rid of empty rows   */
  246.      {
  247.       if ( ! emptyrow( i, nsteps ) )              /* all entries in row = 0  */
  248.       {
  249.        for( j=i+1; j<nmetab; j++ )
  250.         if ( emptyrow( j, nsteps ) )              /* j = first non-empty row */
  251.         {
  252.          rowsw( i, j, nsteps );                   /* switch r